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Abstract Interactions between swimming bacteria have led to remarkable ex¬ 
perimentally observable macroscopic properties such as the reduction of the 
effective viscosity, enhanced mixing, and diffusion. In this work, we study an 
individual based model for a suspension of interacting point dipoles repre¬ 
senting bacteria in order to gain greater insight into the physical mechanisms 
responsible for the drastic reduction in the effective viscosity. In particular, 
asymptotic analysis is carried out on the corresponding kinetic equation gov¬ 
erning the distribution of bacteria orientations. This allows one to derive an 
explicit asymptotic formula for the effective viscosity of the bacterial suspen¬ 
sion in the limit of bacterium non-sphericity. The results show good qualitative 
agreement with numerical simulations and previous experimental observations. 
Finally, we justify our approach by proving existence, uniqueness, and regu¬ 
larity properties for this kinetic PDF model. 

Keywords effective viscosity • kinetic models • bacterial suspension • 
asymptotic analysis 


Contents 


1 Introduction . 2 

2 Model for Semidilute Bacterial Suspensions. 4 


M. Potomkin 

Department of Mathematics, The Pennsylvania State University, University Park, PA 16802 
E-mail: mup20@math.psu.edu 

S. D. Ryan 

Department of Mathematical Sciences and Liquid Crystal Institute, Kent State University, 

Kent, OH 44240 

E-mail: sryanl8@kent.edu 

L. Berlyand 

Department of Mathematics, The Pennsylvania State University, University Park, PA 16802 
E-mail: berlyand@math.psu.edu 







2 


Mykhailo Potomkin et al. 


2.1 Definition of the effective viscosity for a suspension of point force 


dipoles. 7 

3 Conditions imposed to derive an explicit formula for the effective viscosity 9 

3.1 Separation of variables. 9 

3.2 Existence of a steady state Pd(d) . 9 

3.3 Px(x) is constant in the 2 :-direction. 10 

4 Derivation of asymptotic expression for small B . 10 

4.1 Evaluation of Fourier transforms . 12 

4.2 The form of asymptotic expansion in P. 13 

4.3 Contribution at 0{B) . 14 

4.4 Contribution at O(P^) 15 

5 Explicit formula for the effective viscosity . 17 

5.1 Mechanisms required for the decrease in the effective viscosity ... 18 

5.2 Effective noise conjecture. 18 

5.3 Comparison to numerical simulations . 19 

6 Global solvability of the kinetic equation. 21 

7 Conclusions. 28 

A Appendix: Explicit form of integral terms Ij from (39) . 30 

B Appendix: Justification of the representation formula (17) . 33 


1 Introduction 

Bacterial suspensions exhibit remarkable macroscopic properties due to the 
emergence of self-organization among its components. In particular, interest¬ 
ing effective properties such as enhanced diffusivity, the formation of sustained 
whorls and jets, and the ability to extract useful work among other results have 
been recently observed for suspensions of bacteria, such as Bacillus subtilis [40, 
37,22,34,8]. The striking experimental observations on the effective viscosity 
provide the motivation for studying a suspension’s effective properties; namely, 
the observation of a seven-fold reduction in the effective viscosity of a suspen¬ 
sion of swimming B. subtilis [35] . This reduction is observed below 2% volume 
fraction typically referred to as the dilute regime where bacteria are far apart 
and essentially interact with the background fluid only. With the assumption 
of no interbacterial interactions, this regime has been studied analytically in 
recent works (e.g., [32,17,15,16]). There bacterial tumbling was introduced 
in order for the formula to predict a decrease in the effective viscosity [16]. 
However, in the absence of tumbling (e.g., for anaerobic bacteria) the decrease 
is still observed experimentally [35]. It was shown recently in [29] that in¬ 
terbacterial interactions substantially contribute to effective viscosity and an 
estimate for this contribution was given. Rigorous analysis of this contribution 
and its corresponding effect on the effective viscosity of the suspension is the 
main component of this paper. 

We begin with an individual based model (IBM) previously introduced 
in [29,28], which has been successfully used to capture the decrease in the 
effective viscosity and other collective phenomena. Such suspensions, where 
interbacterial interactions play an important role and are modeled as a sum 
of pairwise interactions, are referred to as semi-dilute. Our goal is to identify 
the underlying mechanisms that contribute to the decrease of the effective 



















Effective Rheological Properties in Semidilute Bacterial Suspensions 


3 


viscosity in this concentration regime. The main tool we employ is a kinetic 
theory derived from this individual based model. 

The purpose for employing a kinetic approach is to replace a large sys¬ 
tem of coupled differential equations by a single continuum partial differential 
equation with respect to a probability distribution of bacteria positions and 
orientations. Note that it is natural to consider probabilistic quantities since 
the main focus of this work is the study of the effective properties. The main 
computational advantage of the kinetic approach is that the number of bacte¬ 
ria N does not increase the complexity of the problem [39,5]. Namely, the PDE 
could be solved numerically with a fixed spatial or temporal grid independent 
of N. In addition to the ability to consider many different initial conditions 
at once, another advantage to introducing this probabilistic framework is to 
consider the limiting regime as —>■ oo, the so-called mean field limit. More 
information on kinetic equations can be found in the seminal works of the 
1970’s [24,6,10] or more contemporary reviews [7,19,25,9]. 

Significant difficulty in the analysis come from the incorporation of inter¬ 
actions. First, they appear in the kinetic equation as a non-local term due 
to the fact that the suspension of interacting bacteria is generally described 
analytically by configurations of all bacteria. Second, the main interactions 
that are taken into account are hydrodynamic, which diverge as bacteria ap¬ 
proach one another as the square of their distance. This results in a singular 
kernel in this non-local term. Thus, the kinetic equation consists of a nonlocal, 
nonlinear PDE due to the presence of interactions. 

Using a kinetic approach, the main result of this paper is an explicit asymp¬ 
totic formula for the effective viscosity with interbacterial interactions taken 
into account. The formula reveals the physical mechanisms necessary for the 
decrease in effective viscosity observed experimentally. To achieve this result 
we first find the steady state solution of the kinetic equation and then use this 
solution to compute the effective viscosity. For completeness, we also establish 
the well-posedness of the kinetic equation. 

This paper is organized as follows. Section 2 begins by introducing the 
individual based model under consideration for a semi-dilute bacterial suspen¬ 
sion. From this, the kinetic equation for the orientation distribution is formally 
derived. The reason we begin with the IBM is that the effective properties of 
a suspension are derived from knowledge of microscopic configurations, which 
is transferred from the IBM to the kinetic model. In Section 3 we introduce 
the main conditions under which we derive the asymptotic formula for the ef¬ 
fective viscosity and discuss their physical significance. Section 4 contains the 
derivation of the asymptotic steady state solution to the kinetic equation for 
the orientation distribution in the limit of small non-sphericity. The effective 
viscosity from the asymptotic formula is then compared to the same quantity 
computed from direct simulations of the individual based model in Section 5. 
The important physical mechanisms for the decrease in viscosity are identified 
and the orientation distribution is compared to the results of previous works 
in the dilute case. In addition, the normal stress differences and relaxation 
time are considered. The existence, uniqueness and regularity properties of a 
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solution to the kinetic PDE are proven in Section 6. Finally, we formulate our 
conclusions and outline potential future investigations in Section 7. 

2 Model for Semidilute Bacterial Suspensions 

We begin by introducing the coupled PDE/ODE system governing the fluid 
and bacteria dynamics respectively. Each bacterium is represented as a point 
force dipole. One force represents the bacterium’s propulsion mechanism (e.g., 
flagellar motion) and the other is the opposing viscous drag exerted by the 
bacterium’s body on the fluid. This approximation has been experimentally 
verified by observing the flow due to a bacterium (e.g., Bacillus subtilis) in a 
fluid and comparing it to that of a force dipole [11]. 

As a bacterium swims through the fluid its trajectory may be altered 
through interactions with other bacteria and the background flow. At every 
moment in time a bacterium propels itself in the direction in which it is ori¬ 
ented. If one bacterium comes into close contact with another, then a collision 
can occur altering the bacterium’s position. This is modeled by an excluded 
volume potential. Finally, the flow itself has an impact on a bacterium tra¬ 
jectory through the ambient background flow and the sum of flows induced 
from the propulsion of all the other bacteria on its position. To make these 
ideas more concrete we now introduce an individual based model (IBM), which 
governs a bacterium’s position and orientation. 

We consider N bacteria with the position of the center of mass of the 
Ah bacterium x® = (ccbybz*) and orientation d® = {d\,d\^d\). A bacterium’s 
translational velocity is derived from a balance of forces due to self-propulsion, 
collisions, and the flow field acting on the position of the bacterium. A bac¬ 
terium’s orientation velocity is derived from a balance of torques in the form 
of Jeffery’s equation for an ellipsoid in a linear flow with additional terms 
due to the flows generated by the other bacteria in the suspension [20]. Thus, 
the equations of motion for bacterial positions x and orientations d originally 
introduced from first principles in [29] are 



( 1 ) 





where Vq is an individual bacterium’s swimming speed and B is the Brether- 
ton constant which takes into account the geometry of the bacterium’s body 
(B <C 1: near spherical, i? « 1: needle-like). The externally-imposed planar 
shear flow contributes to each bacterium’s motion through the fluid velocity. 
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= (0,7a:, 0)^, as well as its effect on a bacterium’s orientation through the 
vorticity Wq ® = Vx x and rate of strain Eq ® ^ (VxU^*^ + (VxU^*^)^). 
Here TV is a white noise and we let I? ~ be the diffusion coefficient. This 
order of D will be used throughout this work and represents the idea that the 
random motion present in the system has a greater effect the more elongated 
a particle is. 

The additional terms in Jeffrey’s equation (2) beyond the contribution from 
the background flow are due to the vorticity and rate of strain generated 
by the j-th dipole on position of the i-th dipole 

u>i = Vx X E^' = i (VxU^ + (VxU^y) . 

Each of these terms depends on the fluid velocity u-^ , which is governed by 
Stokes equation and will be described in greater detail below. 

Remark 1 The equations of motion (l)-(2) are a 5V coupled system of or¬ 
dinary differential equations in comparison to the dilute case studied in [16] 
where there were only two ODEs governing the evolution of a single bacterium 
in an infinite medium (only depending on a single bacterium’s orientation). 
Thus, the semi-dilute system of equations adds a greater complexity than the 
dilute case previously studied. 


The use of Stokes equation to model the fluid is justified by estimating the 
Reynold’s number. Based on the typical size £0 ~ 1 and swimming speed 
Vq ^ 20 fim/s of a bacterium, in addition to the typical dynamic viscosity 
r]o ~ 10“^ Pa • s and density p ^ 10^ kg/m^ of the suspending fluid, the flow 
has a Reynolds number Re around 2 x 10“^ ^ 1. Thus, inertial effects can 
be neglected. Also, it is assumed that a steady-state flow is established on a 
timescale much smaller than the characteristic timescale, which is the time for 
a bacterium to swim its length. 

The flow at the position of bacterium i due to bacterium j is given by 
u^ (x*, d-1) = u(x'’ — x*, d-l) where u(x, d) is a solution of the Stokes problem 

( ?7oAxu(x, d) - Vxp(x, d) = Vx • [D(d)(5(x)], x G 

< Vx ■ u(x, d) = 0, X G (3) 

[ u(x, d) —)• 0, jxj —>• 00, 


where 770 is the ambient fluid viscosity and p is the pressure. The dipole tensor 
D = {Dim} is given by 


Dim{d) ■= Uo 




( 4 ) 


where Uq is the strength of the dipole referred to as the dipole moment. For 
pushers, bacteria that propel themselves from behind such as B. subtilis, Uq < 
0. Equation (3) has an explicit solution: 


Ufe(x, d) 


1 

871770 


3 3 

EE Dlm{(i)Gkl,m{^), 

1—1 m —1 


( 5 ) 
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where 0fci(x) = (^ + is the Oseen tensor. 

Remark 2 In order to study the role of interactions in semi-dilute suspensions 
it is natural to deal with a point representation of swimmers such that the 
whole suspension is modeled by points interacting in the fluid. In our paper, 
a swimmer is represented by a point force dipole with the dipole tensor (4). 
In general, for a given model of a swimmer, such a point representation can 
by found as the second order term in the multipole expansion, see [21]. We 
note that all results of this paper such as the asymptotic formula for orien¬ 
tation distribution and effective viscosity can be easily modified to semi-dilute 
suspensions with swimmers whose dipole tensor is different from (4) . 

In order to analyze the system (l)-(2), the associated kinetic theory for 
the probability density of bacterial configurations (positions and orientations 
of each bacterium) is studied. In general, to derive the corresponding kinetic 
equation one assumes that initial conditions are random. Then each sum in the 
equations of motion is a sum of identically distributed random variables. The 
key step in the formal derivation of the kinetic equation is replacing all sums 
in the equations of motion by their expectations [26,39,18]. This allows one 
to replace all the sums representing interactions by integrals with respect to a 
probability density function Pit, x, d) of finding a given bacterium at position 
X with orientation d. 

By replacing the sums with integrals in the system (l)-(2) and enforcing 
conservation of probability, a standard Fokker-Planck equation describing the 
evolution of the density P is obtained 

dtP + V^-(yP) + Vd-{f2P)-DAdP = 0, ( 6 ) 

where the translational and orientation fluxes are defined by 

V(x,d) := Vbd-h [ uP(x',d')dx'dS'd'+ u®‘^(x), (7) 

I CL I JVl 

I2(x, d) / (a; + BE, P(x', d'))x'd^d' + a;^G(d) -k PE^^(d). (8) 

I Cl I Js^ 

Here < •, • > denotes the duality with respect to the L^-norm, Vl '■= [—L, L]^, 
and we neglect the Lennard-Jones term F due to the fact that collisions only 
play a small role at low concentrations. The functions it, o’, and E under the 
integral sign depend on x — x', d and d', and they are defined as follows 

■= ■ [(fid - //3)^l(x)], 

u;(x,d,d') := -id X [Vx X u(x,d')] , (9) 

E(x, d, d') := —d x [d x i7x(u(x, d'))d], 

where Px(u) := g(VxU-|- [Vxu]^) represents the symmetric gradient and I is 
the identity matrix. Also, and E^®(d) are defined in the same way 

as (9), but with the fluid velocity u replaced with the background flow 
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Remark 3 Since a;,E ~ the integrals with respect to the spatial vari¬ 

ables must be considered in the distributional or principal value sense (which 
are equivalent here). Namely, 



where 



The orientation vector d € 5^ can be represented by two independent 
angles in spherical coordinates 


d := (cos a sin/?, sin a sin/3, cos/3) = {di,d 2 ,dz)^ 


( 10 ) 


for azimuthal angle a G [0,27r) and polar angle /3 S [0, tt) with unit basis vec¬ 


tors a := (—sina,cosa, 0) and jd := (cos a cos/3, sin a cos/3, — sin ^) respec¬ 


tively. Here one must be careful to note that the divergence and the Laplacian 
in orientations (the Laplace-Beltrami operator) in (6) are taken over the unit 
sphere. In particular, for any field A = A(d) the following definition holds 



( 11 ) 


where = A ■ a, Ap = A -13, and Vd is the classical gradient. 

2.1 Definition of the effective viscosity for a suspension of point 
force dipoles 

To define the effective viscosity consider the contributions to stress: (i) due to 
dipolar hydrodynamic interactions 



l,m = 1, 2, 3, 


depending only on each particle’s orientation [3] and (ii) due to soft collisions 
(the excluded volume constraints) 



depending only on the relative positions of each bacterium [41]. Both are 
combined to form the total stress due to interactions first used in [29,28]. 
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We assume that all bacteria are in the volume Vl at any instant of time. The 
bacterial configurations are denoted by x := (x^,x'^) and d := (d^,d^). 

The ultimate goal is to compute the effective viscosity due to hydrody¬ 
namic interactions at low concentrations for comparison with experimental 
observation [35] and numerical simulations. At lower concentrations cj), where 
the striking experimental decrease in the effective viscosity was observed, the 
contribution due to collisions is relatively small and for the proceeding analysis 
will be neglected 

i:(x, d) = E‘^{d) + « A'^(d), for cj) small. (12) 

The exact concentration interval where the formula (12) works well will be 
determined later by comparison with direct numerical simulations of the sus¬ 
pension. 

Thus, it is sufficient to restrict attention to the density of orientations 
denoted fd(d) defined as 

Pd{d) := ^ [ P(x, d)dx, where [ Pd(d) = 1. (13) 

JVl 

For comparison with experiment, the main quantity of interest is the shear 
viscosity or component 771212 of the fourth order viscosity tensor relating the 
stress to the strain, henceforth, denoted as r). We define the effective viscosity 
as the averaged ratio of the corresponding components of the stress and strain 
tensors 

[ —^(x,d)dxd5d = ^ / Ei^{d)Pd{d)dSd, (14) 

Vo \yL\JvLJs^ 7 iJs^ 

as in [29,28]. Here p = N/\Vl\ is the mean concentration or number density. 

The following nonlinear, nonlocal integro-differential equation describes the 
evolution of the orientation density Pd{t, d) 

dtPd{t,d) = -Vd • (< >x -Pd(i, d)), (15) 

where < >x= ^ Jy^flPx{t,x)dx, Ct contains the background flow and 

interaction terms 

n{t, X, d) = -t [ [ (a;-kPE,P(t,x',d'))dx'dd'. 

A|Kz,| 752 

Equation (15) is obtained by integrating (6) in x and dividing by N. 

Remark 4 In this work, lower concentrations of bacteria are considered where 
the primary contribution to the effective viscosity from interactions is the dipo¬ 
lar component of the stress, which only depends on the set of bacterium 
orientations. Thus, the x equation will not factor into the final formula; how¬ 
ever, F is the force associated to a truncated Lennard-Jones type potential 
imposing excluded volume constraints. For more information on its definition 
and why it is needed for global solvability see [28] . This quantity still remains 
in the original coupled ODE system used for simulations to ensure that parti¬ 
cles remain a finite distance apart avoiding an artificial divergence in the fluid 
velocity u ^ l/|x® — x-^ | (see section 5.3). 
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3 Conditions imposed to derive an explicit formula for the effective 
viscosity 

To calculate the effective viscosity we impose three conditions to make the 
system more amenable to mathematical analysis. 


3.1 Separation of variables 

In this paper only small concentrations are considered where collisions are 
not important, yet the flow of each bacterium affects all others. The bacteria 
are at large distances apart and, thus, since the background flow provides the 
major contribution to bacterial motion, then distributions of positions and ori¬ 
entations become essentially independent of one another. This can be justified 
from the experimental work of Aranson et al. (e.g., see [37,2]). Henceforth, it 
is assumed that the positions and orientations are decoupled. 

Condition (Cl): The density P(x, d) can be written as 

P(x, d) = Px(x)Pd(d) (separation of variables), (16) 

where Pd(d) = P(^, d)dx and fg2 Pd(d)dSd = 1. Here N is the number 

of bacteria, supp(Px(x)) C Vl, where the spatial density Px(x) can be found 
by Px(x) = Jg2 P(x,d)dS'd. 

This condition is used twice. First, the effective viscosity at low concentra¬ 
tion only depends on the orientation (see Remark 4). Thus, using condition 
(Cl) an explicit equation for the evolution of the orientation distribution can 
be derived from (15). Second, V formally contains diverging integrals (e.g., 
J J FdxdSd since F ^ |x|“^^), which will no longer be present in the equation 
for the orientation distribution Pd(d) allowing for further mathematical anal¬ 
ysis. It will be observed at the end of this work that the asymptotic expansion 
for Pd(d) depends on Px(x) through the coefficients, thus all the information 
about spatial patterns is preserved. 


3.2 Existence of a steady state Pd(d) 

A steady state solution to (15) is defined as follows 
Definition 1 Pd(d) is called a steady state solution to (15) if it solves 
0 = —Vd • ^< D >x Ai(d)^ . 

To compute time independent effective viscosity we impose the following 
condition. 

Condition (C2): There exists a nontrivial steady state solution to (15). 
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First, note that there is no trivial steady state unless B = 0 in which case 
we find the uniform orientation distribution (d) = ^. This can be obtained 
both in the limit as S —?► 0 in the asymptotic results derived herein for Fd(d) 
and from observing that the trivial steady state would be a constant satisfying 
the constraint fg 2 Pd{d)dSd = 1- One still needs to prove the existence of a 
steady state in the general case B ^ 0. The condition (C2) can be formulated 
as a theorem and its proof may be the topic of a future work. Here we remain 
focused on the study of the effective viscosity. 


3.3 T’x(x) is constant in the z-direction 

We assume that Px(x) is constant in z for the case of the planar shear back¬ 
ground flow under consideration in this work. This is consistent with past 
numerical observations by Ryan et al. [29] and experimental observation in 
[38] since the suspension remains below any critical concentration for three- 
dimensional collective motion. Also, collective motion even in full 3D experi¬ 
ments and simulations in planar shear flow has been observed to be essentially 
2D in the shearing plane [38]. Thus, following experimental observation, we 
assume the same. 

Condition (C3): The density T’x(x) is constant in z. 

The condition (C3) essentially follows from the physical setup of the quasi- 
2D thin film suspension. In Appendix B we show that the condition (C3) leads 
to the following representation formula for the Fourier transform of the spatial 
distribution F’[Px]: 

iF[P^]f = S{h)PUki.k2). (17) 

Here k = {ki, k 2 , k^) is the Fourier variable, and ^12(^1, ^2) is a smooth func¬ 
tion defined in k-space independent of k^. 


4 Derivation of asymptotic expression for Pd for small B 

In this section, an expression for the orientation distribution Pd(d) is derived. 
Since (15) is a nonlinear integro-differential equation it is challenging, in gen¬ 
eral, to find an analytical solution. Thus, we look for Pd(d) by asymptotic 
expansion in the limit of small non-sphericity (P <C 1). This will allow us to 
apply analytical techniques and derive an expression, which will provide phys¬ 
ical insight into the mechanisms contributing to the decrease in the effective 
viscosity. 

Rewrite the equation for the orientation density Pd(d) (15) as (the argu¬ 
ment t is suppressed for simpler notation) 

dtPd + Vd • [(o;®^ -b BE^^)Pd] + I Vd • (/>P(x, d))dx = 0, (18) 
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where 

/?(x, d):=j^ f {uj + BE,(19) 

I til J52 

Herein fi will denote the component of the orientational flux 17 due to inter¬ 
actions. Observe that the u> and E are functions of x — x', d, and d'. 

Using Condition (Cl) defined in (16) we obtain a closed form equation for 
a steady state Bd(d) (provided that P^ is given): 

0 = Vd- [(u;^^ + BE^G)Pd(d)] 

+ X d')Px(x)Pd(d)) dBd-dx. (20) 

The first term in (20) is the contribution due to the background planar shear 
flow: 

Vd ■ [(w^‘^(d) -h BE^'^(d))Pd(d)] = -^^sin^/3sin2Q:Pd(d) 

+ ^(l + Pcos2a){5„Pd(d)} (21) 

'7 B 

+ sin2asin2/3{cl,3Pd(d)}. 

The second term in (20) is the contribution of hydrodynamic interactions be¬ 
tween bacteria. Notice the convolution form of the nonlocal terms in the spatial 
variable. In the next section, the Fourier transform will be utilized to compute 
quantities necessary to derive the formula for the effective viscosity. Specifi¬ 
cally, using tools such as Parseval’s Theorem, one can take the spatial integrals 
and consider them in Fourier space where they will prove easier to analyze. 
After using the separation of variables (16), the density will be expressed in 
terms of the Fourier frequencies k. 

The main goal for the remainder of this section is to write the system in 
a convenient form for using the Fourier transform. This idea follows naturally 
from the aforementioned observation that all the interactions terms take the 
form of a convolution. Introduce the Fourier transform C(k) := P[Px](k): 

= (22) 

Define H(x—x', d, d') := u;(x—x', d, d')-|-PE(x—x', d, d'), then the following 
equalities hold 

< H*Px,Px >x=< P[H*Px],P[Px] >k=< P[H],(P[Px])' >k, (23) 

where * and P stand for convolution and Fourier transform, respectively. The 
first equality is Parseval’s identity and the second is the fact that the Fourier 
transform of a convolution is the product of Fourier transforms. Thus, one can 
rewrite equation (20) in the following form 

Vd - [(a;®^ + PE^^)Pd(d)] 

+ / Vd • {Pd(d)Pd(d') < P[H](P[Px]) 2 >k} dPd' = 0. (24) 
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In order to compute -F[H] one must first understand how the Fourier transform 
acts on the fluid velocity u and its derivatives. 


4.1 Evaluation of Fourier transforms 

In order to analyze (24), an analytical expression for the Fourier transform 
F[H] = E[aj] + i3F[E] is needed. Both terms depend on the fluid velocity u 
dehned by (3). Recall the dipolar stress 

r(x, d) = D(d)(5(x) = Uo{dd - //3)5(x). (25) 

Then the Stokes equation in (3) can be written as 

- 77o4\xU + VxP = Vx • T'(x,d), Vx-u = 0. (26) 

Denote the Fourier transform of a function f{x) as 

f{k) = F[f] (k) = I e-*(‘^-)/(x)dx, 

and compute the Fourier transform of u and the symmetric gradient Ilx(u). 
Proposition 1 Let u be a solution of (3) and let S be defined by (25). Then 
(*) r(d') = Uo (d'd'* - //3), 

(in) E[i:)x(u)] = - ^ flkl^iikk* -2kk*i:kk* + Ikpkk*^) . (28) 

Here * denotes the transpose. 

Proof The part {i) follows from the fact that the Fourier transform of 6- 
function is 1. 

We split the proof of {ii) into two steps: First, we find the Fourier transform 
of the pressure p, then by using the first equation in (3) we find u. 

Step 1: Evaluation ofp = F'[p]. By taking the divergence of (26) in x we obtain 

Z\xp = Vx-(Vx-i7). (29) 


Observe that 

F [A^p] = -|k|2p(k), F [Vx • (Vx • r)] = J S : V^e-*‘^'"dx = -r(k) : kk*. 

Substituting these formulas into (29) we obtain — |kpp(k) = —T’(k) : kk*, 
and, thus, we find an expression for the Fourier transform of the pressure p: 

p(k) = ^^(k) : kk*. 


(30) 
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Step 2: Evaluation of u = r'[u]. Return to Stokes equation (26) and observe 
that 

rjoF [Z\xu] = -ryo|kpu(k), F [Vxp] = ikp(k), 

F [Vx • r] = ir(k)k. 

Using these relations one finds that ? 7 o|kpu(k) + *kp(k) = ii7(k)k. After 
rearranging the terms and using (30) we complete the proof of (ii). 

To prove (iii) we first observe that F [_Dx(u)] = | (uk* + ku*). Plug the 
Fourier transform of u from (ii) into this expression to find 

F[D^(u)] = ^ (uk* +ku*) 

1 / kk* ~ ~ kk* \ 

= - W - lijp - lijp)) ^ 

Use the fact that E is symmetric (E = E*) to complete the proof of (iii). 

Remark 5 It is easily seen that F [_Dx(u)] does not depend on |k|, since F [_Dx(u)] 
can be rewritten as 

r-rr. / M 1 f, k k* 2 k k* ~ k k* ^ k k* ~ 

F[nx(u)j - 

This subsection is concluded by summarizing the analytical expressions for 
the two main components of Up] = U [a;] + BF [Ej: 

U [E] = -d X (d X U [U>x(u)] d) = F [i:>x(u)] d - dd*F [i:>x(u)] d (31) 

U [w] = -^d X F [Vx X u] = -^d X [-ik x F[u]], (32) 

where F[u] and F[Fx(u)] are given by Proposition 1. 


4.2 The form of asymptotic expansion in B 

Recall the steady-state Liouville equation (24) with the background terms 
substituted in: 

0 = sin^ (3 sin2aFd(d) -|- ^(1 -|- F cos 2a)daFd(d) 

'7 B 

+ sin 2a sin 2/35^Fd(d) 

+ ■ {^d(d)Fd(d') < F[H], (F[Fx]) 2 >k} dS^r (33) 

We consider the asymptotic expansion in the Bretherton constant, F ^ 1, for 
the orientation distribution, Fd(d), up to the second order: 

Fd(a, /?) = F^°^ (a, /3) + F^'^ (a, /3)F -f F^'^ (a, P)B^ + 0(B^). (34) 







14 


Mykhailo Potomkin et al. 


Substituting (34) into (33) we get different equations at different orders of 
B. It is straightforward that P) = ^ (surface area of the unit sphere 

is 47r) solves the equation at order 0(1). We want to consider the asymptotic 
expansion about the uniform distribution because it has been extensively doc¬ 
umented in theory and experiment that as the bacterium bodies become or 
spherical {B —>■ 0), then the distribution in angles is uniform [29,16]. In the 
next two subsections, the linear order term P^\a,j3) and quadratic order 
term P^ ’{a,j3) are computed. 


4.3 Contribution at 0[B) 


First, notice that Vd • a;(x — x', d, d') = 0. Indeed, this follows from (11) since 
a; • d = 0 and the classical divergence of a; with respect to d is zero (note 
that a; = d X A, where A — Vx x u does not depend on d). This observation 
implies Vd • .F[H] = BVd ■ .^[E]. 

Using this equality and expanding the divergence under the integral sign 
we rewrite (33) as follows: 

0=1 [i?sin(2a) sin/3 {cos^dpPd — 3sin/?Pd) + (1 + Pcos(2a)) d^Pd] 

+ Pd(d')Pd(d)(Vd • (P[E(d)])(P[Px])^)kd5d' (35) 

+ Vd[Pd(d)]Pd(d')(P[H(d)](P[Px])2)kd5d-. 

The first integral at 0{B) is 

. (F|E(d)|)(F|PJ)>),<iS,„ (36) 


By switching the order of integration and noting fg 2 SdSd' = fg 2 bfo[d'(d')* — 
I/‘S\dSd' = 0 we obtain that (36) is zero using (31) and (28). 

Since both Vd[Pd(d)] and PE are of the order 0(P), the second integral 
in (35) at 0(P) is VdPd^^(d)(P[a;](P[Px])^)kdS'd' which is also 

zero due to fg 2 Uo(d'd' — I/i)dSd' — 0. 

Thus, the integral terms do not contribute to equation (35) at order 0{B), 
and it has the following form: 



-3P 


( 0 ) 


sin(2Q!) sin^ /3 -|- daPd^ 


(37) 


After substituting pj^^ 


^ and solving (37), one finds that 




3 

— sin^/3cos(2a). 
Stt 


( 38 ) 
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Since the integral terms are zeros at order 0{B), the contribution due to in¬ 
teractions does not appear at order 0{B) and thus the only contribution is due 
to the background flow. 

It will be shown later that up to 0{B) the contribution to the effective vis¬ 
cosity by the bacteria is zero. This will shed light on the fact that interactions 
are necessary to see the decrease in the effective viscosity and the background 
flow alone is insufficient. Note that even though this is the contribution due 
to the background flow the strain rate 7 is not present. Therefore, the magni¬ 
tude of the flow will not have an effect on the longtime limit of the effective 
viscosity at 0{B). However, once the terms at the next order are computed 
one observes a competition develop between the background flow and the flow 
due to inter-bacterial interactions. In this case the magnitude of the shear 7 
becomes important. 

4.4 Contribution at 0{B^) 

Consider terms in (35) of order 0{B‘^)-. 


0 = ^ sin(2Q;) sin f cos f dp Pj^\d) -^ si 

+ l9aPd\d) + ^ cos(2a)ac,id^^(d) 


Y sin(2a) sin^(/3)P^^^(d) 



(39) 


Denote the four integral terms in equation (39) by H, I 2 , I 3 and I 4 , respec¬ 
tively. The following equalities hold: 


40^77oA^|yL| 

I2 = I3 = 0, 

, 3C/o 


(H sin^ (3 cos(2a) -I- C sin^ (3 sin(2Q;)) 


10^77oIV|I4| 


I?sin(2a) sin^ /3, 


where constants A, C, and D are defined as follows 


A -.= ^ fsin^(20)P^2^^dkdO, C := J sin{4e)P?2k^dkd0, 


D := J cos(d) sm{6)Pi2k^dkd9. 


( 40 ) 
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Here P 12 is from (17), and we use spherical coordinates in the Fourier space 
{k = |k|, 6, (j>). The calculations of h can be found in Appendix A. 

After substitution of the expressions for each b, we get the following equa¬ 
tion for 


0 = ^ sin(2a) sin/?cos^5/3Pj^^(d) — ^ sin(2Q;) sin^(/3)Pj^^ (d) 

+ 1 cos(2a)aaPy^(d) 


(41) 




C/o 


m^T^oN\VL\ 

3Po 


(A sin^ /3 cos(2a) -|- C sin^ j3 sin(2a)) 
Psin^ /?sin(2Q;). 


io^inoN\VL[ 

Based on the form of the equation (41), the following representation is used 
to hnd P^\d): 

P^^^(d) = Cl sin"^/3cos(4a)-|-(72 sin^/3cos(2a)-I-Ca sin^/3sin(2a). (42) 


In order to find each Ci substitute (42) into (41): 


0 = 


£ - 


- 7 C 2 + 


7 C 3 + 


sin(4a) sin^ ■ 

UqC 

AOtttjoNIVlI 10^77oA^|Vl|J 




407r7;oA(|Fi| 


sin^ P cos(2q;) 


WqD 


sin^ P sin(2a). 


Since the factors are linearly independent, each coefficient is zero and, thus, 
we find the C,’s: 


Cl = 


167r ’ 


^ ^ Uo{C+l2D) 

^ 4077rr/oAr| Vy I ■ 


C3 = - 


UqA 

4077r77oAf|y£, | ' 


Using these coefficients one obtains an explicit formula for the orientation 
distribution up to 0{B^): 


1 3 

PAa,P) = -sin^/3 cos(2q;)P-I- 

47r Stt 


IGtt 


sin"^/3cos(4a) 


~ Po ^ I P cos(2q;) 


4077r?7oA(|y£,| 

UoA 2 


(43) 


4077r77oA(|UL| 


sim P sin(2Q;) 


-hO(P^). 


Formula (43) is the main result of Section 4. Since A, C, and D contain P 12 , 
all the spatial information is embedded in these coefficients. In particular, 
we found the lowest order (in B) contribution of hydrodynamic interactions 
to the Pd(d) occurs at 0{B^). In the following section, the contribution of 
hydrodynamic interactions to the effective viscosity is computed as well as the 
change in the effective normal stress coefficients. The combination of these two 
quantities will describe the total effect of hydrodynamic interactions on the 
rheological behavior of the bacterial suspension. 
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5 Explicit formula for the effective viscosity 

Using the expression for the orientation distribution, rd(d) defined in (43), 
and the formula for the effective viscosity for dipoles in a suspension (14), we 
compute the contribution to the effective viscosity due to interactions: 


^int V-Vo 
Vo 


U^B^P^A 

75^^Trr]Q 


(44) 


where A = -^A ^ 0 ( 1 ) and the equality holds up to order 0{BA- The 
quantity 77 '"* behaves like in concentration (cf. [4] where an expansion 
for the effective viscosity to order two in concentration is derived for pas¬ 
sive spheres corresponding to pairwise interactions). As an additional check of 
consistency, consider the dimensions of the final quantity. The dipole moment 
[h^o] = ; both the Bretherton constant B and A are dimensionless, the 

concentration/number density [p] = 7 ^, the ambient viscosity [ryo] = 7 ^, and 
the strain rate [ 7 ] = ^ resulting in ry™* being dimensionless. 

In addition, the orientation distribution Pd(d) from (43) can be used 
to compute the effective first and second dipolar normal stress coefficients 
Ni 2 = and N 23 = to investigate the effect of hydrody¬ 

namic interactions. The main advantage of the mathematical model is that 
the computation of the effective normal stress coefficients is straightforward in 
contrast to experiment where its measurement can be quite complicated [13]. 
These coefficients can provide important information about the suspension. 
For example, the ratio of the first normal stress to the viscosity determines 
the effective relaxation time [13]. Also, phenomena such as extrudate swelling 
[1] and secondary flow [27] are important in many technological applications. 
A simple calculation shows that 


Ni 2 = 
Ni3 = 


y'd Y'd 
^11 ~ ^22 


y<d Y'd 

^22 ^33 


UqP 

UoP 


2 2Uop{C+12D) ^^ 

5 JbyTrryo 

1 ^ Uop{C+12D) ^^- 
5 TSyTrryo 


(45) 

(46) 


The approximations are valid for B ^ 1, so for pushers {Uq < 0) N 12 > 0 
and N 23 < 0 where as for pullers {Uq > 0) N 12 < 0 and N 23 > 0. Both 
results are consistent with the predictions in [16,31] while providing additional 
information about the concentration dependence. The effective normal stress 
coefficients grow linearly with concentration in the presence of interacting 
bacteria; however, the fact that the normal stresses of active suspensions are 
non-zero in the case of a planar shear flow indicate the emergence of non- 
Newtonian behavior. One sees in (45)-(46) that as the shear rate 7 —> 00 the 
normal stresses approach zero indicating the dominance of the background 
flow on the suspension overwhelming any contribution from interactions. 
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5.1 Mechanisms required for the decrease in the effective viscosity 

In this subsection, the mechanisms that lead to a decrease in the effective vis¬ 
cosity are investigated. These same mechanisms are shown in [30] to be respon¬ 
sible for collective motion and large scale structure formation in suspensions of 
pushers. Our mathematical analysis provides insight beyond experiment. For¬ 
mula (44) reveals that elongation of bacteria, self-propulsion, and interactions 
are all required to observe a decrease in the effective viscosity; namely, for 
spherical bacteria {B = 0) the net change in the effective viscosity is zero. In 
addition, active bacteria are required, since f7o ~ /p = 0 results in no change 
in the effective viscosity where fp is the propulsion force. Finally, if the spatial 
density T’x(x) is near uniform, then A = f sm^( 20 )Pi 2 dk « 0 resulting in 
no change in the effective viscosity. 

In the limit 7 —> oo the contribution to motion of bacteria due to shear 
dominates the contribution due to interactions with fd(d) maximized at a = 
7 r /2 and /3 = 7 r /2 (alignment with y-axis). This is analogous to the passive 
case where bacteria in a planar shear flow tend to align with the direction 
where the fluid exerts the least amount of torque on the bacterium body. 
Therefore, confirming our main conclusion that in order to exhibit a decrease 
in the effective viscosity active, elongated bacteria whose interactions result in 
a non-uniform distribution in space are needed. 


5.2 Effective noise conjecture 

In this subsection, the results herein involving a semi-dilute suspension of 
point force dipoles are compared to the previous result for a dilute suspension 
of prolate spheroids with propulsion modeled as a point force [16]. Thus, the 
only contribution to bacterial motion is the background flow. In [16], finite 
size bacteria are taken as spheroids with a point force (6 function) accounting 
for self-propulsion. In addition, each bacterium experiences a random reorien¬ 
tation referred to as tumbling. Biologically tumbling corresponds to a reori¬ 
entation of a bacterium in hopes of finding a more favorable (nutrient rich) 
environment. Typically in experiment this is observed when the concentration 
of oxygen is low. Thus, bacteria enter a more dormant state resulting in a 
lower swimming speed and an increased tumbling rate [36] . 

Since only the term containing A contributes to the effective viscosity, one 
can choose to match the coefficient of this term 



sin^/3 sm(2a) + O(B^) 


with the corresponding coefficient in the derivation by Haines et al. [16], which 
is quadratic in the diffusion strength D. To make the formulas for the effective 
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viscosity identical, the strength of the effective noise/diffusion (tumbling) is 
chosen to be 


— 15r?o7^ + \/225r?Q7'^ — A^B'^^'^p'^Uq 

D ■= - -X---- > 0, 

UABpUo 

(since Uq < 0 for pushers). Observe that D, chosen in this way, depends only 
on the physical parameters present in the problem and the same effective 
viscosity as the dilute case studied in [16] is found. This D is referred to 
as the effective noise and the phenomenon where stochasticity arises from a 
completely deterministic system is called self-induced noise. A future work may 
seek to explain this phenomenon rigorously using mathematical analysis. One 
heuristic idea is that the periodic (deterministic) Jeffrey orbits are destroyed 
by interactions resulting in stochastic behavior. 

Some conclusions about this effective noise can be made that ensure its 
consistency with physical reality. As bacteria become spheres B ^ 0, D ^ 0 
resulting in no change in the effective viscosity consistent with [16]. Also as the 
strain/shear rate 7 —>■ oo, I? —>■ 0. This is physically intuitive, because as the 
shear rate becomes large its contribution dominates that due to hydrodynamic 
interactions resulting in behavior that resembles that of a passive suspension. 
Thus, the contribution to the effective viscosity due to hydrodynamic interac¬ 
tions is zero. Finally, we compare our results with direct simulations for the 
coupled PDE/ODE system composed of Stokes PDE (3) and (l)-(2). 


5.3 Comparison to numerical simulations 

In this section, the accuracy of the derived formula is tested by comparing it 
to recent numerical simulations. The numerical procedure is outlined in [29]. 
These simulations are parallel in nature allowing them to be carried out on 
GPUs for greater efficiency. 

Figure 1 shows how both the formula and numerical computations of vis¬ 
cosity change with bacterium shape as all other system parameters remain 
fixed. Here shape is accounted for through the Bretherton constant B = 557^ 
where b is the length of the major axis and a is the length of the minor axis of 
the ellipsoid representing a bacterium. First, notice that in both the formula 
and numerics the contribution to the effective viscosity due to hydrodynamic 
interactions decreases with B (increasing in magnitude). This is due to the 
fact that as bacteria become more asymmetrical as i? —>■ 1 the inter-bacterial 
hydrodynamic interactions have a greater effect on alignment. This alignment 
increases the magnitude of the dipolar stress leading to an even bigger decrease 
in the effective viscosity. The agreement between the analytical formula and 
numerical simulations breaks down as B becomes large, but this is expected 
due to the fact that the asymptotic formula is valid in the regime where B 1 
(small non-sphericity). 
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0 0.1 0.2 0.3 0.4 

Bretherton Constant. B 

Fig. 1 Comparison of the formula for the effective viscosity with numerical simulations as 
bacterium shape changes through the Bretherton constant B for a fixed volume fraction 
'P = .02 and shear rate 7 = .1. The vertical bars represent the error in the numerical 
approximation. Error in the analytical solutions comes from the numerical estimation of A. 


Figure 2 shows how both the formula and numerical computations of vis¬ 
cosity change with the concentration of the suspension as all other system 
parameters remain fixed. It is seen that as concentration increases the effec¬ 
tive viscosity decreases. This can easily be explained by the fact that as the 
concentration increases, the motion of bacteria begins to be dominated by 
inter-bacterial hydrodynamic interactions. This leads to collective motion of 
the bacteria in the suspension, which subsequently decreases the viscosity. The 
two results begin to diverge near volume fraction <P « .02. The reason the nu¬ 
merical simulations do not decrease as much is that collisions are taken into 
account. It was shown in [29] that the stress due to collisions is a positive con¬ 
tribution to the effective viscosity that is not captured by the formula. This 
contribution begins to become important beyond the dilute regime > 2%). 



Fig. 2 Comparison of the formula for the effective viscosity with numerical simulations as 
the volume fraction ^ changes for a fixed shape B = .2 and shear rate 7 = .1. 












Effective Rheological Properties in Semidilute Bacterial Suspensions 


21 


Figure 3 shows how both the formula and numerical computations of vis¬ 
cosity change with the shear rate of the background flow in the suspension as 
all other system parameters remain fixed. As expected when the shear rate 
is large in both the analytical formula and simulations, the decrease in vis¬ 
cosity due to hydrodynamic interactions is negligible. This is due to the fact 
that the background flow dominates motion of bacteria wiping out the ef¬ 
fects of inter-bacterial interactions and stopping any collective structures from 
forming. When the shear rate is too small the effective viscosity becomes un¬ 
bounded. This makes sense given that at small shear rate the system becomes 
almost non-dissipative and thus the effective viscosity is not well-defined. This 
can easily be seen by noting that the viscosity is the ratio of the stress over the 
strain and when the strain is essentially zero the effective viscosity becomes 
unbounded. All three plots show good qualitative agreement with each other, 
experimental observation, and physical intuition. 



Fig. 3 Comparison of the formula for the effective viscosity with numerical simulations as 
the shear rate 7 changes for a fixed volume fraction ^ = .02 and shape B = .2. 


6 Global solvability of the kinetic equation 


In this section, we study solvability of the main nonlinear integro-differential 
equation (15) governing the evolution of the orientation distribution. Primarily 
we are interested in existence, uniqueness, and the regularity properties of 
solutions of (15). 

First, we note that (15) is an equation of the form: 


dtPd = -Vd 


Us^ 


K{d, d')Pd{d')dSd' +kid) 


Pd I + D^dPd- 


(47) 


Indeed, one can obtain (15) by substituting 

A:(d, d') = a;(d, d') -f BE(d, d'), k(d) = a;^'^(d) -1- BE®^(d). 


( 48 ) 
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Both K and k from (48) are infinitely smooth functions of d. Therefore, in 
this section we consider (47) for the general case of smooth K and k. 

We follow the standard procedure for the analysis of the well-posedness 
of the evolution PDEs (e.g., see [12,23,14]). In particular, we introduce the 
notion of a weak solution. By (s S K) we denote the corresponding Sobolev 
spaces. 

Definition 2 For T > 0, the function / which belongs to space T-L given by 
H = L^{{0,T),HHs‘^)) n (49) 

is a weak solution of (47) if for almost all t G [0,T] and all h G 


(dt/,h) = -D(Vd/,Vdh) + (/, 


K{d,d')fdSd' +Hd) 


■ Vdh), 


(50) 


where (•, •) is the duality product for distributions on the unit sphere S^. 

Remark 6 According to the well-known embedding (see [33]) the fact that a 
weak solution / belongs to % implies that it is continuous with respect to 
t G [0, T] with values in L^(5^), i.e., / S C'([0, T]; L^(5^)). 

Definition 3 A function / G C{[Q,T]; L^(S^)) is called positive in distribu¬ 
tional sense if 

{f,h)>0 (51) 

for all t G [0, T] and all h G C{S^) such that h(d) > 0 for all d G 5^. 

The following theorem is the main result of this section. 

Theorem 1 Assume fo G K G (7^(5^ x 5^), k G (7^(5^) and T > 0. 

Assume also that /o is positive in the distributional sense. Then the following 
statements hold: 

(i) There exists the unique weak solution of (47) / on interval [0,r] such 
that f\t=o = fo- The weak solution f is positive. It continuously depends 
on initial conditions, i.e., there exists a positive constant (7 > 0 such that 

sup 11/(1) _ /(2)||^2(52 ) < (7||/,^^) - /,$^^||l 2(52), (52) 

tG[0,T] 

where /(i) and /(^) are weak solutions with initial conditions /(i) |t=o = 
and /(^)|t=o = fo^\ respectively. 

(a) For all s > 0 if fo G then f G C([0, T]; i7®(iS^)). 

Iffo G C-(52), thenf G Ci[0,T];C°°{S^)). 
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(Hi) For all s > 0 if fo G then for all m > 0 and t > 0: 

||/(t)||^.+r.(52) < C ^ , 


(53) 


where the constant C depends only on ||/o||ff»(52)) s, and m. In particular, 


f GC{{o,^y,HP{s^)) 


for all p G Z. 

Proof 

STEP 0. {Preliminaries) Consider spaces of functions with mean zero: 

L^S^) := L^{S^) n {/ : (/, 1) = 0} H^{S^) := H%S^) D {/ : (/, 1) = 0} . 

Note that for / G L^{S^) 


if, 1)= [ fdS^. 

Js^ 

We use ||Vd/||L2(52) as a norm in H^{S‘^). 

In this proof we assume that f ^2 fo^Sd = 1. Consider the “mean zero” 
component of the solution /; namely, g := / — ^. If / is the weak solution of 
(47), then g satisfies 


d 

dt 


{9,h) 


-n(Vdg, Vd/i) + {^+g, f K{d, d')g{d')dSd' ■ Vd/i) 

47r Jg2 




/52 


X(d,d')d5'd' +k(d) 


Vd/i) 


(54) 


for all h G H^{S^). Existence, uniqueness, and continuous dependence on 
initial conditions will be proven for g, which is equivalent to the proof of the 
same properties for /. 

Below C denotes a positive constant and it may change from line to line. 

STEP 1. {Local existence) Let Ejv be the space spanned by the first N eigen¬ 
values of the Laplace-Beltrami operator Z\d, and let iljv be the orthogonal 
projector on the space E^- Introduce the Galerkin approximation g^, which 
is the solution of the following equation: 


d 

dt 


{9^,h) 


-D{Wd9^, Vd/i) + {^+9^', f K{d, d')g^{dydSd' ■ ^dh) 

dvr J^2 




U52 


K{d, d')dSd' + k(d) 


Vdli), 


(55) 


for all h G Em, and g’^\t=o = IlNgo, where go := fo - 

In a standard manner, the problem (55) can be interpreted as a system of 
N ODEs, and its solution g^ exists for t G [0,tjv) for some tM > 0. Taking 
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h = in (55), using the Cauchy inequality, and the boundedness of K and 
k we obtain 

^Il5'^lli2(52) + ^ C* ^1 + 115^11^2(52)^ . (56) 

In the inequality (56) the constant C does not depend on N. This implies that 
exists for 0 < t < tQ where to may be chosen independently from N, and 

llff'^(t)||L2(52) <C, 0<t<to, (57) 

The bound (57) gives that the RHS of (56) is estimated by a constant inde¬ 
pendent from N. Then by integrating (56) in t we get 

r\\9^\\lHsndt<C. (58) 

Jo 

Take h G L^(0, to; in (55), integrate in t, and use the Cauchy inequal¬ 

ity, {u,v) < C'||M||iji(52)||?;||^-i(52), and the Minkovsky inequality to obtain 


pto 


{dtg , h)dt < C 


-I 1/2 


|^i( 52 )dt 


Therefore, 

r \\dtg^fH-^^s-)dt<C. (59) 

^0 

From bounds (57), (58), (59) and the following relation which holds for all h 
from % 




{dtg, h)=- {g, dth)dt + {g{to), h{to)) - (g(0),h(0)), 


we obtain that there exists g G H such that (up to a subsequence) 


5^-5 mL°-{0,to-,L\S^))nL\0,to;H\S^)), (60) 

dtg^^g in L^iO, to; (61) 


In particular, weak convergences in (60) and (61) imply strong convergence in 
C{[0,to];L^S^)). Thus, 

5|t=o = lim 5^|t=o = lim iIjv5o = 5o, 

N—>-oo N—¥oo 


and 



K{d,d')g^id')dSd' 


K{d,d')g{d')dSd' in Ci[0, to]; L^S^)). (62) 


To complete the proof of local existence we need to show that g solves (54). 
To this end, consider (55) with h = w{t)ho, where ho G Em, M < N and w{t) 
is arbitrary smooth function of one argument t. Integrate this equation in t 
over the interval (0, to), and pass to the limit N ^ oo {M is fixed) using (60), 
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(61) and (62). Since w{t) is arbitrary we obtain that (54) is satisfied for all /ig 
from the space UmEm which is dense in Therefore, g solves (54) for 

all h G 

Thus, we constructed a function g that is a weak solution of (54) defined 
on the time interval 0 < t < tg. 

STEP 2. {Uniqueness & continuous dependence on initial conditions) 
Consider g^^'^ and g^‘^\ weak solutions of (54) defined on the time interval 
[0,to] with initial data and gl^\ respectively. For both i = 1 and i = 2 if 
one substitutes h = into the equation (54) written for g^'‘\ one obtains by 
using the same arguments as for (57) that 

E (63) 

i^l,2 

where the constant C depends on initial data gj^'^ and the parameter D only. 

By subtracting equation (54) written for from equation (54) written 
for g^^^ we get the following equality 


{dtu,h) = -£)(VdU, Vdh) + 


/S2 


K{d, d')udSd> 


,Vdh) 


Us^ 


is^ 


K{d,d')dSd' +1^ 

K{d,d')g^^'>dSd' 




is^ 


K{d,d')udSd' 


Vdh) 

Vdh) 

Vdh). 


(64) 


By taking h = u, using the Cauchy inequality, and (63) we obtain 

^Iklli2(52) + ^ll^llffl(52) < C'||u|||2(52). 

This inequality implies that \\u{t)\\‘j^ 2 (^g 2 -^ < e^*||^t(0)|||2(52), and, thus, 


- go^h^s^)- (65) 

Again, the constant C depends on initial data g^'^ and the parameter D only. 

The inequality (65) implies that a weak solution of (54) continuously de¬ 
pends on the initial data. In particular, uniqueness holds: if gj^^ = g^\ then 
from (65) it follows that the corresponding solutions g^E amt gU) coincide. 

STEP 3. {Regularity of weak solutions) 

Consider a weak solution g and assume go G H^{S^) that s G Z+. Such a 
weak solution exists due to STEP 1, and it can be approximated by Galerkin 
approximations g^ which follows from uniqueness proved in STEP 2. 
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By substituting h = {—Ad)^g^ into the equation (55), using the Cauchy 
inequality and (57) we obtain 

~^\\ 9 ^\\‘ h ={ s ^) + ^Il5^llff»+i(52) ^ C'(ll5^llff»(S2) + 1); (66) 

where the constant C depends on | 1 ( 7 oIIl 2 ( 52 ), || 5 ollff=( 52 ) and the parameter D. 
In the same manner as for (57), (58) and (59) it follows from ( 66 ) that 

is bounded in n 

dtg^ is bounded in (0, to; (S^)). 

Hence, g S L^{Q,to]H‘^"^^{S^))r\L°°{QAo\H‘^{S^))r\H^{Q,to-,H’^~^{S^)). The 
standard embedding theorem (e.g., from [33]) implies g G C{[0,to]; 

STEP 4. {Positivity of weak solutions) 

Consider / = ^ + 5 , a weak solution of (47). Assume first /o G H'^{S‘^) 
and /o(d) > 0. Then / belongs to (^([O, to]; C'^(‘5^)), and thus / is a classical 
solution of (47): 


dtf = DAdf - F • Vd/ - (Vd • F)/, 

where F(d) := Jg 2 K{d,d')f{d')dSd' + k(d) G (^([O, to]; ^^(5^)). Consider 
/ := where uj := max jV^ • Fj. Then / solves the following equation 

[0,to]x52 

dtf = DAdf - F ■ Vd/ + (w - Vd • F)/. 

Since w — Vd • F > 0 the weak maximum principle for parabolic equations 
applies for /, and, thus, / > 0. 

Consider the case of /o G Lf{S^), which is positive in the distributional 
sense. Then we can approximate /o by positive fo G H'^{S‘^) in the space 
F^(5^). Denote by f^ solutions of (47) with initial data f^. Then by (65) we 
can pass to the limit A^ —>■ 00 in the inequality 

{f^{t),h)>0 

for all 0 < t < to and h G C{S^). Thus, the function /, which is the solution 
of (47) with initial data /o, is positive at least in the distributional sense. 

STEP 5. [Global existence) 

Consider /o = ^ + ffo G F^(5^), which is positive in the distributional sense. 
Functions / and g are weak solutions of (47) and (54), respectively. We want 
to prove in this step that the time interval on which / and g are defined can 
be extended from [0,to] to [0,T] for any given T > 0. 

First, observe that 


/ f{t)dSd = [ 


fodSd = 1 . 


52 


52 
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From the equality above and positivity of / established in STEP 4 we obtain 

II/(0IIli(52) = 1- 

In particular, since |g| < |/| + ^ we have 

[ K{d,d')g{d')dS^, < C(||/(t)|Ui( 52 ) + 1) = 2C. (67) 

Js^ 

Substitute h = g into (54), use the Cauchy inequality and (67) to obtain 

^Il5lli2(52) + D\\g\\jji(^g2'^ < C (^||5||i2(52) + . 

Then the L^-norm of the weak solution is bounded on all bounded time inter¬ 
vals [0,T]: 

j^max^ 115(^)11^2(52) < + 1). 

Thus, global existence follows. 

STEP 6. {Instantaneous regularity) 

Consider positive /o S i7®(5^) and the corresponding weak solution / = 
^ + g oi (47). According to STEP 3 / G L^([0, T]; i7®+^(5^)) and, thus, 
/ G iJ®+^(5^) for almost all t > 0. Hence, there exists t > 0 arbitrarily 
close to 0 such that f{i) G i7®+^(5^). Then by uniqueness and STEP 3, 
/ G C([t,T];i7®+^(iS^)). We can choose i arbitrarily small and T arbitrarily 
large (due to global existence proved in STEP 5). By repeating the same 
arguments for s -I- 1, s -I- 2, and so on we get 

f GC{0,+^;HP{S^)) 


for all p G Z. 

Next we prove (53) by induction with respect to m. Substitute h = (—Z\(i)^(7-|- 
<(— Z\d)®^^5 for t > 0 in (54) and use the Cauchy inequality to obtain 

^ ^Il5lljj»(52) + y^ll5lljj» + l( 52 )^ 

+ Y ^II6 'IIa=+i( 52) + Y^ll5ll^s+2(52)^ < C'll5o||ffq52)(l -I- t). 

Using the Poincare inequality Hfflijjs+fc(52) < ||5lliji.+fc+i(52) we obtain 

A'’(52 ) + ■^^I!5II^s+i(s2) < C'||po|Ih=(52)(1 -I- 1). 

Thus, the base of induction is shown 


11511^0 + 1(52) < C'||(7o||^o(52) 


1-H 


1 

t 


( 68 ) 
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Finally, to get the inequality (53) at the order m + 1 we use the inequality 
(53) at order m between times t/2 and t and (68) between times 0 and t/2: 

|| 5 (t)||^=+™+i( 52 ) < C\\g llff»+i( 52 ) ^ 

< C'll5o||ff=(52) (^1 + ■ 

Thus, (53) is proved by induction. 

STEP 7. {Proof of Theorem 1) 

(i) Existence of a weak solution of (47) for arbitrary T > 0 is proved in 
STEP 5. Uniqueness is proved in STEP 2. To prove continuous dependence 
on initial data on arbitrary time interval [0, T] one needs to repeat all 
arguments in STEP 2 replacing tg by T. Positivity is proved in STEP 4. 

(ii) This part is proved in STEP 3, if one replaces tg by T. 

{Hi) This part is proved in STEP 6. 

□ 


7 Conclusions 

In this paper, the derivation of a formula for the effective viscosity formally 
derived in [29] was made rigorous and an additional term in the asymptotic ex¬ 
pansion for the effective viscosity was derived (now up to 0{B^)). This formula 
revealed the physical mechanisms responsible for the decrease in the effective 
viscosity confirming the prior formal calculation. Namely, hydrodynamic in¬ 
teractions, an elongated body, and self-propulsion are required to observe a 
decrease. These features are all present in the bacteria Bacillius subtillis used 
in the experiments of Aranson et al. [37,34,35,38,36], which motivated this 
study of the effective viscosity. In addition, an interesting phenomenon was 
uncovered: the emergence of self-induced noise where a completely determin¬ 
istic system governed by interactions resembles a random system for certain 
regimes of the physical parameters. The explicit analytical formula for the 
effective viscosity derived herein showed good qualitative agreement with sim¬ 
ulations and experiment. This paper also establishes the global solvability of 
solutions to the PDE kinetic equation governing the evolution of the bacterium 
orientation density. In order to derive the formula for the effective viscosity, 
the existence of a steady state was assumed and then computed asymptoti¬ 
cally. Rigorously proving the convergence to a steady state distribution may 
be the subject of future work. 
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A Appendix: Explicit form of integral terms li from (39) 


We will need the following technical Lemma: 


Lemma 1 Assume that A is a 3 X 3-matrix that is independent of the orientation vector 
d. Then 

Vd • [d X (d X yld)] = 3(d, AA) - TvA. (69) 


In particular, if 


A = 


A C 0 
C -AO 
0 0 0 


(70) 


then 


Vd ■ [d X (d X Md)] = A sin^ ^ cos(2o) -h C sin^ /3 sin(2cK). 


(71) 
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Remark 7 Recall that Vd denotes the spherical gradient in orientation d, and Vd denotes 
the classical gradient in vector d (e.g., see (11)). 


Proof Using the well-known vector identity a x {b x c) = 6(a, c) — c(a, b) and the relation 
(11) we obtain 


Vd ■ [d X d X ^d] = Vd ■ [d(d, ^d) - ^d] 

= Vd ■ [d(d,^d) - Ad] - ^ ||d|5(d,^d) 
9|d| I 


(72) 


|d|3(d,yld)} 


ldl=l 


Here d = d/|d|. The orientation d is a unit vector, but in order to relate the classical and 
the spherical divergence we need to calculate the derivative in |d| at |d| = 1; thus, consider 
d different from unit magnitude. Also, note that d does not depend on |d|. 

One can easily verify that 

Vd ■ [d(d, ^d) - ^d] = 3(d, ^d) + d ■ V(d, ^d) - Tr(^) 

= 3(d, ^d) + d ■ 2Ad - Tr{A) 

= 5(d,^d)-Tr(A). (73) 


and 


■—{|d|5(d,yld)-|d|3(d,^d)} 


|d|=i 


= -2{d,^d). 


(74) 


Substituting (73) and (74) into (72) we obtain (69). The formula (71) follows directly from 
substituting (70) into (69). 

□ 


Next we evaluate integral terms R introduced in Subsection 4.4. 
The integral term R. This integral is defined by 

and due to (28) and (31) can be written as: 


II 


-1 


SjjttWIVlI 


d X d X / MF[P^fdk \ d 


Here k := k/|k|, and the 3x3 matrix A4 is defined by 


M := J'kk* - 2kk*J='kk* + kk* 


where 



31/o 

Stt 


- Stt 
15 


0 

0 


0 

Stt 

15 

0 


0 

0 

0 


% 

5 


10 0 
0-10 
0 0 0 


(75) 


Substituting T into the expression for fA one finds that fA equals to 

— 2k^ + 2k^k2 —2kfk2 + 2k\k2 k\k^ — 2k\k^ + 2k\k2kz 

— 2k^k2 + 2kik^ ~‘^^2 ~ ‘^^ 1^2 + ^^2 —^ 2^3 — 2kfk2k3 + 2k^k3 

kik^ — 2kfk^ + 2kik2k^ —k2k^ — 2k^k2k^ + 2fc2^3 —2kfk^ + 2k2k^ 


where /ci, ^ 2 , ks are components of k. 
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Next we use the condition (C3) from Subsection 3.3 written in terms of the representa¬ 
tion formula (17) to obtain 


j MiF[P^]fdk = j >!U3^oA2(|k|fci,|k|fc2)|k|2d|k|de, 

where 


•^ 1 * 3=0 — 


■ 2 fc 2 [l_fc 2 _,_fc 2 ] 

— 2 A:ifc 2 [A;i — fc|] 
0 


— 2 fcifc 2 [fci — ^ 

— 2 fc|[l + fcj — fcj] 0 
0 0 


4^3 A:| —2fcifc2[fci — ^ 

— 2 A:iA: 2 [fcj — fc|] 0 

0 0 0 


sin2(20) -isin(40)O 
-isin(40) -sin2(20) 0 
0 0 0 


Here variables fci and k 2 are expressed in polar coordinates fci = cos(fl) and k 2 = sin( 6 ). 

Note that the matrix in the equality above is of the form (70) and, thus, we may apply 
(71): 

= 77;— I 0 cos(2a) + C sin^ 0 sin(2o)) . 

407r77o-/V| Vx| 

This is the desired formula for Ii. 


The integral term I 2 . We need to obtain that I 2 = 0. This holds provided that 


J{F[u;]{F[P^]f)^dSa, = 0 . (76) 

52 


The integral in the RHS of (76) by using inverse Fourier transform can be written as 

— -dx J J 'X u(x — x^, d^)d5'(j/j> 


The integral in curly braces is zero due to 

f u(x,dOdd' = [ f/o(d'(dO* -I/3)dSa' 

Js^ Us^ 

Thus, (76) holds and I 2 = 0. 


: Vx^ = 0 . 


(77) 


The integral term I 3 . To prove that I 3 = 0 we can use the same arguments as for I 2 . Indeed, 
I 3 vanishes provided that 


|(F[E](F[PJ)2)kd5d- =0. 


52 

One can easily verify this equality by using the inverse Fourier transform and the identity 
(77). 

The integral term I 4 . This integral can be written as 

[ pi'^(d')(PM(F[Px])")krfSd.. 

N\Vl\ 752 


According to (32), the formula for F[a;] is the following 
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Recall that k = k/|k| = (/ci, /c 2 , fca). Using (75) we obtain 


M-.= [ 
Js^ 



lOr? 


X 


[k X kk*) Pic] 


k2ks 
kiks 
— 2fci/C2 


In the same manner as we analyzed Ii, we use the condition (C3) from Subsection 3.3 
written in terms of the representation formula (17), the form of orientation d in spherical 
angles (10), and polar angle 6 for ki = cos 9 and k 2 = sin0: 

l4 = ■|MU 3 ^oA" 2 (|k|fci,|k|fc 2 )|k|"d|k|de, (78) 

where 

— sin a sin 0 
cos a sin 0 
0 

Next we find Vd[Pd^^(d)]. Using (38) and the definition of the spherical gradient Vd^ 


A^U 3 =o = ^ sin( 2 e) 


VdP = 


- °in|^j Qc»P + cos(o) cos(i3)dpP 
s°n(g) + sin(a) cos(l3)dpP 
— s\yi{0)dpP 


we obtain that 


VdP^^d) 


sin OL sin(2a) sin 0 + cos cx. cos(2q) sin 0 cos^ 0 
— cos Oi sin(2Q;) sin 0 + sin a cos(2q:) sin 0 cos^ 0 
— sin^ 0 cos 0 cos(2q:) 


Substituting this equality into (78) one obtains the desired formula for I 4 : 

3Uo 


I 4 = 


-Dsin( 2 Q!) s\n^ 0. 


107r?7oN|VL| 

This concludes the evaluation of integral terms U for f = 1,4. 


B Appendix: Justification of the representation formula (17) 


Consider the spatial distribution P:x.(x^y,z) = CLx{^)Pi 2 (_^-,y)i where 



{-L.L), 

{-L.L), 


(79) 


and we choose cl = ‘ily/irL. The distribution Px satisfies the condition (C3), i.e., its support 
does not depend on 2 :. 

Our main goal of this subsection is to obtain a representation for the Fourier transform 
ofPx: 

P[Px]=/ x(2)e“3"dfc3A2(fcl,fc2) = -^"^'^^^A2(fcl,fc2). (80) 

J-L VttL kz 

For an arbitrary continuous function (f) the following convergence holds: 


1 /■ / sin(fc3L) 

Al J V A 

-L 


2 

(f>{k3)dk3 —>■ (pio). 


(81) 
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From (80) and (81) it follows that for large L we have 

{F[P^]f ^S{k3)Ph{kuk2), (82) 

which justifies (17). Note that due to our choice of cl it follows from fy^ = N and 

TV ~ that Pi 2 y/Jj. 

It is also interesting to calculate the coefficient A defined by (40) for the spatial distri¬ 
bution Px(x) = ^x{^)xiy)xi^) which is uniform in Vl- Then 

=vz - 5ikuk2)L^/\ (83) 

In this case the coefficient A is of the order It is responsible for the decrease in 

viscosity. Namely, for fixed number density p = N/L^, the Bretherton constant B, the dipole 
moment Uq and the strength of the background flow 7 it follows from (44) that ~ Ajh^. 
Then due to (83) 

^ 0 as L ^ 00 . (84) 

Therefore, A = AL~'^I'^ can serve as a measure of the deviation of the spatial density Px(x) 
from uniform. 





